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Abstract 



Recent results, extending the Schmidt decomposition theorem to wave- 
functions of pairs of identical particles, are reviewed. They are used to 
give a definition of reduced density operators in the case of two identical 
particles. Next, a method is discussed to calculate time averaged entan- 
glement. It is applied to a pair of identical electrons in an otherwise empty 
band of the Hubbard model, and to a pair of bosons in the Bose-Hubbard 
model with infinite range hopping. The effect of degeneracy of the spec- 
trum of the Hamiltonian on the average entanglement is emphasised. 

1 Introduction 

Schmidt decomposition 

Assume that the wavefunction ?/>(xi,X2) describes two distinguishable particles. 
Then there exist orthonormal bases of wavefunctions m (xi) and x m (x2) and 
coefficients p rn > such that ip can be written as a single sum 



This result is known as the Schmidt decomposition theorem. See for instance 
PQ, Theorem 2.7. The reduced density matrices for each of the particles are 
then given by 




(1) 
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Indeed, one verifies that for any one-particle operator A 



TrerA = }^p m {(p m \A\(p m ) 



rn 



* jan.naudts@ua.ac.be 
ttobias.verhulst@ua.ac.be 



1 



= "y]p m {(l>m ® Xm\A ® l\(/> m <8> Xm) 
m 

= (4) 

and similarly 

TttA = O A|V>- (5) 

The knowledge of the coefficients p m suffices to calculate the von Neumann 
entropies 

£(ip) = — Trainer = — Trrlnr = — ^~~^p m lnp m . (6) 

m 

The latter quantity is a measure for the entanglement of the two particles. 
Identical particles 

Recently El Hj , the previous result was generalised to pairs of identical par- 
ticles, described by a wavefunction ip in a Fock space. Let &t(0) and b((j>) be 
the creation and annihilation operators for a particle with wavefunction </>(x). 
Let |0) denote the vacuum state. Then for each two-particle wavefunction ip in 
a Fock space there exists an orthonormal basis of wavefunctions <p m (x) in the 
one-particle Hilbert space and coefficients p m > such that 

i/> - -^=Y,VP^bH&n)bH<j> m )\0), (bosons) (7) 

mions). (8) 

m 

If the dimension of the one-particle Hilbert space is odd then the latter expres- 
sion does not involve all of the basis vectors <j) m . 

The physical interpretation of this result, in the case of bosons, is that with 
probability p m the two particles are both in the same state with wavefunction 
4> m - In the fermionic case, one of the particles is in the state (j>2 m , the other in 
the state 4>2m+i- It is then obvious to define reduced density matrices a and r 

by 

d = T = y^p m \<t>m) {<t>m\ (bosons), (9) 

m 

O = 2^P2m|^2m) (<t>2m\ and T = 2 P2m+1 1 4>2m+l ) (02m+l|, 
m m 

(fermions) (10) 

By convention, P2n+i = P2n in the latter case. 

In the fermion case these density matrices are far from unique since for any 
pair (/>2m,fen+i the two basis vectors may be interchanged. Nevertheless, the 
resulting values of the von Neumann entropies of a and r are always the same. 
Hence, in all cases the quantity 

£(lp) = - ^p„lnp n (11) 
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can be used as a measure of entanglement. 

In the next sections we reproduce the proofs of (0 \TU\i and show that the 
eigenvalues p n of the reduced density matrices can be calculated without actually 
performing the generalised Schmidt decomposition. In this way the quantifica- 
tion of the entanglement of a pair of identical particles is more easy than in the 
case of distinguishable particles. 

Linear entropy 

Even the simplified method to obtain the eigenvalues p n may be too difficult 
for analytical treatment. For this reason we will make use of the linear entropy 
instead of the von Neumann entropy (llip . It is still a measure of entanglement 
[5], and is given by 

£iM=X>»(1-1>tO (12) 

n n 

For similar reasons the von Neumann entropy has been replaced by the linear 
entropy in other papers as well, for instance in O [3 [8] . 

The simplification arises as follows. Let p be a density matrix with eigen- 
values p n . Then it is often feasible to calculate p 2 by matrix multiplication 
while the calculation of phip usually requires diagonalisation of p. Also calcu- 
lating the trace of p 2 is usually a feasible task. The linear entropy £i(ip) is then 
obtained as 1 — Tr p 2 . 

Average entanglement 

A final simplification comes from averaging the linear entanglement. In prin- 
ciple, the entanglement of two particles depends on time. Rapid fluctuations 
of entanglement have been reported to occur in vibrational modes of triatomic 
molecules [5] , and between electrons of Rydberg molecules [TU] • They have been 
studied in theoretical models such as the Dicke model [TT] , a model of coupled 
kicked tops [12] , the Harper Hamiltonian [13] , a dimer model [14] , Bose-Einstein 
condensates [TS] . Hence it is obvious to study the time average of the entangle- 
ment. In |16j it is shown how to replace the time average of non-linear quantities 
such as the entanglement by ensemble averages. This was applied by the present 
authors to study the entanglement of distinguishable particles [T7] . 

Overview of the paper 

The next section recalls known results about symmetric and anti-symmetric 
matrices. Proofs are given in the Appendix. The theorems of [3] [31 0] are repro- 
duced and the calculation of the average entanglement is explained. Section 3 
discusses the entanglement of a pair of identical electrons in an otherwise empty 
band described by the Hubbard model. Section 4 demonstrates the importance 
of degeneracy of the spectrum of the Hamiltonian for a two-boson model. The 
paper ends with a discussion in Section 5, followed by two Appendices. 
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2 Schmidt decomposition in Fock space 



Known results on symmetric and anti-symmetric matrices 

Remember that a matrix M is normal if it commutes with its hermitean conju- 
gate . The transpose M T of M has matrix elements (M T ) mn = M nm . The 
matrix M is symmetric if M T — M, it is anti-symmetric if M T + M = 0. Any 
matrix with complex entries M can be written as M = V^DU with D diagonal 
and with U and V unitary. This is the singular decomposition of M. A similar 
result for symmetric matrices is the following theorem. It is known as Takagi's 
factorisation theorem - see [TS], or [3J, Theorem 3.4. See [15] . Theorem 5.5.1, 
for the first claim of the theorem. 

Theorem 1 Let be given a square matrix M with complex entries. Then M is 
symmetric if and only if it can be written as 

M = U T DU (13) 

with D diagonal and U arbitrary. The matrix U can be chosen unitary. 

Consider for example the matrix M, given by 

m=(; {). d4) 

It is not normal. Still, there exists a unitary matrix U, namely 

and a diagonal matrix D = [1, 1] + (1 + i)[l, -1]/V5 such that M = U T DU. 
The method to find U is based on the observation that U diagonalises M'M. 
Indeed, one has 

M f M = U ] D^DU. (16) 

This observation is essential for the calculations that follow. 

The analogous result for anti-symmetric matrices is usually formulated for 
matrices with real entries only. For matrices with complex entries it follows 
from Lemma 1 of [2]. As noted in [20] , the Theorem below is known in the 
Physics literature since long — see [21] . 

Theorem 2 Let be given an anti- symmetric matrix M with complex entries. 
Then there exists a unitary matrix U such that M can be written as M = U T DU, 
where D has on each row and each column at most one non-vanishing element. 

If M is anti-symmetric then also D = (U yMW is anti-symmetric. Hence, 
if D has at most one non- vanishing element on each row and each column then it 
can be brought into block-diagonal form with blocks of size at most two, simply 
by swapping the order of rows and of columns. This is, D is a block matrix of 
the form 

D= [Z U Z 2 ,--- ,Z X ,0,0,---], (17) 

with Zj of the form 

*' = ( Z o ) ■ ^ 
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Application to wavefunctions in Fock space 

Take an arbitrary orthonormal basis of wavefunctions to n (x) in a finite dimen- 
sional one-particle Hilbert space. Any two-particle wavefunction ip can be writ- 
ten as 

j> = / , \nnU m ® U n - (19) 

The matrix of coefficients X mn is denoted A. In the boson case A is symmetric, 
in the fermion case it is anti-symmetric. Hence, by the previous theorems there 
exists a unitary matrix U and a matrix D, with at most one non- vanishing 
element on each row and each column, such that A = U T DU. Then one can 
write 

if> — ^2 U rrn D rs U s 

mnrs 

= ^D rs (t) r (g)(t) s , (20) 



with 



b r — U rm ui m . (21) 



Because the matrix D has at most one non- vanishing element on each row and 
each column, the double sum in ([2*0"j) reduces to a single sum. This yields (0 [HI 

num. 

Next observe that 

A* A = {U T DU^U T DU = U ] D^DU. (22) 

Hence, the matrices A^A and D'D have the same eigenvalues. But the eigenval- 
ues of D'D are precisely the coefficients p n appearing in the expression (jTTJ) for 
the entanglement. Hence, in order to calculate the entanglement of two iden- 
tical systems it suffices to expand the wavefunction if) in an arbitrary basis, as 
done in (|19p. Next, the matrix of expansion coefficients A is used to form At A. 
Finally, the eigenvalues p n of the latter matrix are calculated. 

Average entanglement using the linear entropy functional 

If now the linear entropy is used to quantify the entanglement instead of the 
von Neumann entropy then one finds 

e x {il>) = 1- Tr(AtA) 2 . (23) 

Next assume that the basis of eigenvectors tp n diagonalises the Hamiltonian H . 
One can expand an arbitrary wavefunction ip in this basis 

<A = EVp^ x "^' ( 24 ) 

3 

with real phases \n and positive coefficients p n satisfying J2 n P n = With 
each basis vector ip n corresponds an anti-symmetric matrix A^"-* via (|19[) . One 
then obtains 

£i{ip) = 1- J2 VPmPnPrPse l{Xn ~ Xm) e l{x °- Xr) Tr A^ 1 A < n > A W f A « . (25) 
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Assume now that the spectrum of H is non-degenerate. Then the time-average 
entanglement of ip may be calculated as an ensemble average, by integrating 
over the phase factors in the above expression. The result is 

&Jjpj = 1 -^p mPr TrA (m)t A (m) A (r)t A (r) 

- Tr A(m)t A(n) A(n)t A(m) 

+ 5>«> Tr (A (m)t A( m )) 2 
= Si(a) +Si(t) - A, (26) 

a = Jj^AW'aW (27) 

m 

r = J^A^A^ (28) 

n 

A = 1 -J^p^Tr (A< m > f A< m >) . (29) 



with 



Note that A^'aW an d AWaW* have the same eigenvalues. Hence one has 

always Si (a) = Si(t). 

The entanglement £i (ifj) calculated above depends on the choice of the basis 
of eigenfunctions of the Hamiltonian. When the spectrum is non-degenerate 
then these eigenfunctions are unique up to a complex phase factor, which has 
no influence on the entanglement. Hence the problem of non- uniqueness occurs 
only when the spectrum is degenerate. In that case the decomposition (I24[) of 
tjj into eigenfunctions should be replaced by 

i' = J2^'n (30) 

n 

with 

^ = w4l ^ Pn = WnW 2 - (31) 

Here, the F n are the orthogonal projections onto the degenerate eigenspaces of 
the two-particle Hamiltonian. 

Examples of degeneracy are discussed below. 



3 The Hubbard model 

As a first application of our method we consider the average entanglement of 
a pair of identical electrons in an otherwise empty conduction band. A suit- 
able description is given by the one-dimensional Hubbard model. There is an 
extended literature about this model. Its study accelerated after Lieb and Wu 
[2"2ll2"5] showed that its spectrum can be calculated using the Bcthc ansatz. For 
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a review paper see [24 . In our treatment here both electrons have the same 
spin. Hence, the Hamiltonian can be simplified to 

N 

H = -J2 hkfyk, (32) 

where bk is the annihilation operator for an electron at site k and the conjugate 
b\. is the creation operator. The coefficients tjk satisfy 

tjj+i = = 1 and tj_k = otherwise. (33) 

Periodic boundary conditions are assumed, identifying site N with site 0. 

We will show that the average entanglement of the two electrons is a non- 
trivial conserved quantity of this model. 



Entanglement of the eigenvectors 



Consider a wavefunction ip describing two identical electrons, say, both with 
spin up, in an otherwise empty band. Then ip is an eigenvector of H, with 
eigenvalue e, if and only if the anti-symmetric matrix A of coefficients X mn 
satisfies the matrix equation 



TA + AT = -eA. 



(34) 



In the one-dimensional model with nearest neighbour interactions (i.e., t mn = 
A(5 m ^ n +i + <5m+i,n) and with periodic boundary conditions (i.e., tjy-ifl — 
to,N—i — A) the solutions are parameterised with two integers r and s, with 
r^s, and are given by 



f 



Ny/2 



[9{mr + ns) — 9(nr + ms)] 



with 0(m) = exp(2irim/N). The corresponding eigenvalue is then 
E (rs) =-2m{r)-2m{s). 

Note that A( rs ' = -A( sr ). 

With the explicit expression (|35|) it is straightforward to calculate 



^(rs)t^(rs) 



1 



2N 2 



V" [9(tr + ms) - 9[mr + ts)] 



x [9(tr + ns) - 9(nr + ts)] 
1 

2N 



9(ms)9(ns) + 9(mr)9(nr) 



Hence, one obtains 



Tr 



^ (rs)t j^(ra) 



I ^ , ^ l 9{ms)9(ns) + 9(mr)9(nr) 

mn 

x ^9(ns)9(ms) + 9(nr)9(mr) 
1 

2' 



(35) 



(36) 



(37) 



(38) 



7 



One concludes that all two-particle eigenvectors tp^ rs ^ are entangled, 
with £i(V> (rs) ) = 1/2. 

One can do even more. The vectors with components 

u ± = 6{mr) ± 6{ms) (39) 

are eigenvectors of the matrix A^ rs ^A^ s ^ with eigenvalue 1/2. All other eigen- 
vectors have eigenvalue 0. Hence, with the notations of previous sections the 
only non- vanishing eigenvalues are po = Pi = 1/2. The entanglement of the 
two-particle eigenvectors ip( rs \ using the von Neumann entropy, is therefore 

£(V^) = 2(-ilni) = ln2. (40) 



Average entanglement 

Let us now calculate the average entanglement of an arbitrary two-particle wave- 
function. One has 



(41) 



Similarly is 



Tr 



J Y(rs)t^(rs) 



A (rV)t A (rV)l = -L^\^\e( ms )e(ns) + e(mr)e(nr) 

mn 

x \(){ns')9(ms') + 6(nr')6{mr' 
= - [6 SS > + 6 rr > + 6 rs > + 8 sr i] . 



Hence 



S(a) = S(t) = 1 — - ^ PrsPr's' [S aa i + &rr' + 8 rs ' + ^sr'} 
rr' ss' 

Using (|26p and the normalisation condition 

^Prs = 1 



(42) 



(43) 



(44) 



one calculates 



£iW0 = o + 



1 i 

2 + 2 



E^ 



2 ^ ^ PrsPr's' [$ss' 
rr' ss' 

1 X - , 

T, + 2^ PrsPr's' , 



+ S rr i + 5 rs i + 6 sr < 



(45) 



where the summation ' is restricted to the sets of indices rr' ss' satisfying 
r > r' , s > s' , r =/= s, r =/= s' , r' =/= s, r' =/= s' . 
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In the above calculation the degeneracy of the spectrum has been neglected. 
As a consequence, the result is only valid when the projection of ip on any of the 
degenerate subspaces is always parallel to one of the basis vectors ip^^ ■ This 
is not the case in general. The calculation of the entanglement of an arbitrary 
wavefunction is therefore more complicated. We will not treat this general case 
but end this section with an example where degeneracy does not play. The 
complications due to degeneracy will be discussed in the bosonic example of the 
next section. 

Example with N = A 

Take N — 4. This means that the two electrons occupy 4 sites on a ring. The 
eigenvalues are -2,0,2, each twofold degenerate. The corresponding eigenvectors 
are ■0 (1 - 4) and -0 (3,4) , ip (1 ' 3) and ip (2 ' 4 \ and tp (1 ' 2 ^ and V (2,3) - We neglect the 
effect of the degeneracy on the average entanglement with the argument that it 
can be lifted by adding a small perturbation to the model. 
Let 

t/j = VpV> (1,4) + p?A (2 ' 3) - (46) 

Projection of tp onto the eigenspace with eigenvalue -2 gives the former term, 
onto the eigenspace with eigenvalue +2 the latter term. The average linear 
entanglement is 



= l/2+p(l-p). (47) 

4 The bosonic model 

As an example of the bosonic case we consider a model which is similar to the 
boson-Hubbard model [231 [2^1 [27] . 

The bosonic creation and annihilation operators satisfy the commutation 
relations [6j,6j^] = 5j k . The Hamiltonian is given by 

N 

H = +Y, *i* 6 i 6 *- ( 48 ) 

However, unlike in the boson-Hubbard model, the hopping coefficients are not 
restricted to nearest neighbour. They rather satisfy 

t jk = [1-(N- l)e]6 jk + e(l - S jk ). (49) 

This model is known as the Bose-Hubbard model with infinite range hopping 

Degeneracy is very important in this model. Indeed, assume e > 0. Then the 
ground state of the one-particle Hamiltonian is N — 1-fold degenerate. Hence, 
the two-particle system has only three energy levels. We will consider the state 
|1, 1,0,0, ••• ,0), in which the photons are not entangled. Next we calculate 
calculate the average entanglement and show that it tends to 1/2 when the size 
N of the system becomes large. 
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Projection onto invariant subspaces 

The one-particle ground state is TV — 1-fold degenerate with energy 1 — TVe. 
Indeed, one calculates for m ^ n 

H(bl-bi)\0) = Y^(t jm -t jn )b]\0) 

= (l-Ne)(bl-bl)\0). (50) 

TV — 1 of these vectors (b' m — 6^)|0) are linearly independent. The remaining 
eigenstate, orthogonal to the ground states, has eigenvalue 1. Its wavefunction 
is 

1 N 

— J2b]\0) = b\^)\0), (51) 

with 

and ujj the one-particle basis formed by u)j = bj\0). Each of these basis vectors 
can be projected onto this eigenvector 

Wj . = J=0(°)+&. (53) 

The vectors £j are orthogonal to </>(°) and hence belong to the degenerate space 
of eigenvectors. 

The one-particle eigenfunction i/jq determines an eigenstate of the two- 
particle Hamiltonian by 

^(oo) = _L 6 t ( ^(o) )6 t ( ^(o) )|0) . (54) 

The initial state 

i/> = |1, 1,0,0, •••) = b\bl\0) (55) 
is now projected onto the three invariant subspaces by writing it into the form 

|1,1,0,0,---) - l&t(0(°))6t(0(o))| O ) 

+ -^ b t (( ^(o) )6 t (a+6) |o) 

+& t (£i)& t (60|o) 

= V / p^/> ( ° 0) + Vp^tP (11) + Vp^^ {01 \ (56) 

with normalised cigenfunctions i})^^ and normalisation constants p^' T \ It is 
straightforward to find that (see the Appendix B) 

P TV 2 ' P { N + TV 2 ' P TV V TV / ' { ' 
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Entanglement 

Next, one should decompose the eigenfunctions ifj^ 00 \ ip^ 11 ^, ^ 01 ^ into the basis 
vectors 

jk 

The calculation of the matrices A^ 00 \ A^ 11 ), A^ 01 ) is found in the Appendix B. 
— see (|80[ I5T1 [52")l . These are used to calculate the density matrices a r, and 
the average entanglements 

and 

12 8 16 16 
= 2 + N-Nl + Jp-N-i- (60) 

See the Appendix B. The final result is 

*® - 5 + ^M) 2 - (61) 

The average entanglement is always larger than 1/2, is maximal at N = 4 with 
a value of 9/16, and converges as 1/N 2 towards 1/2 for large N. 



5 Discussion 



In a rather long Introduction we have summed up a number of results that 
appeared in the literature. We have reviewed known properties of symmetric 
and anti-symmetric matrices, with proofs in the Appendix A. When applied to 
wavefunctions in a Fock space they lead to the definition of reduced density 
operators for systems consisting of two identical particles. These results are 
known. They generalise the Schmidt decomposition theorem to pairs of identical 
particles. We propose to take this generalised decomposition theorem as the 
basis for defining a measure of entanglement of two identical particles. Up to 
now, many authors have used for identical particles the same expressions as for 
distinguishable particles. This leads to the artificial result that the entanglement 
of two identical fermions is always larger than 1. Subtracting this constant 1 is 
not needed when using the definition (|6|). 

In Section 3, the technique to calculate the time-averaged entanglement is 
explained. The linear entropy is used instead of the von Neumann entropy in 
order to simplify the calculations. The extension of this technique to systems 
of two identical particles is straightforward, using the generalised Schmidt de- 
composition. 

Two applications have been considered, one for fermions, the other for 
bosons. In the Hubbard model the average entanglement of two identical elec- 
trons can be calculated for arbitrary initial conditions. However, in this calcula- 
tion we have neglected the effect of degeneracy of the spectrum of the Hamilto- 
nian. This can be justified with the argument that small perturbations caused 
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by the environment would lift the degeneracy. The average entanglement ob- 
tained in this way is always larger than one half and is a non-trivial conserved 
quantity. In the boson model the degeneracy is much worse, leaving only three 
distinct energy levels. For one particular initial state we have shown that the 
average entanglement can be calculated, taking degeneracy into account. The 
resulting value tends to 1/2 when the size of the system becomes large. 

Related results have been obtained by other authors. Levay et al [5U] consider 
2 fermions in combination with a one-particle Hilbert space of dimension 4. 
Wang and Sanders [21] use the generalised decomposition theorem to decompose 
the state of the system into qubit states. Next they calculate the entanglement 
of one qubit with the others and average over the choice of qubits. Plenio et 
al. [3U1[3TJ[32] have considered the typical entanglement in ensembles of Gaussian 
states. These states differ considerably from the two-particle states considered 
here. Nevertheless, the matrix decomposition theorems might be relevant for 
their context as well. 

Only bipartite entanglement has been considered in the present paper. Mul- 
tipartite entanglement is more complicated and requires additional investiga- 
tion. See for instance [33J GUI GIS] . Neither did we study spatial entanglement 
of identical particles [3J3 [371 [33] , or other measures of entanglement, like con- 
currence [35]. Finally, note that we assume that the time evolution is unitary. 
One expects that, due to interactions with the environment, entanglement will 
fade away. See the review paper [391 . 

Appendix A 

For the sake of completeness, we give here a proof of Theorems 1 and 2. First 
assume normal matrices. 

Proposition 1 // M is normal and symmetric then there exists an orthogonal 
matrix V and a diagonal matrix D such that M — V T DV . 



Proof 

Let be a spectral family in a finite dimensional Hilbert space. Then 

there exists a unitary matrix V and two-by-two disjunct sets I n such that 

E (n) = V ^(n)y (g 2 ) 

where 

= 1 if p = <?e/ n 

= otherwise. (63) 

Note that 

4? =E re /„iVV (64) 

Hence, if E^ is symmetric then all elements E^ 1 are real. This implies that, 
if all E^ n > are symmetric, then V can be chosen orthogonal, i.e. — V T . 
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Let M = J2 n \ n E^ be the spectral decomposition of M with all A„ two by 
two distinct. Then also the E^ are symmetric because of the uniqueness of the 
spectral decomposition and because the transpose of an orthogonal projection 
operator is again an orthogonal projection operator. Hence there exists an 
orthogonal matrix V such that 

M = V T DV with D = 5^A n I<'°. (65) 

n 

□ 

If M is anti-symmetric then 

o = ]Ta„ [e^ + (e^) t 

n 

This docs not imply that the E^ are anti-symmetric (which is impossible for 
a non- vanishing orthogonal projection operator anyway)! Hence a different line 
of reasoning is needed. 

Proposition 2 If M is normal and anti- symmetric then there exists a unitary 
matrix U such that U T MU has on each row and each column at most one non- 
vanishing element. 



(66) 



Proof 

Let M — J2 n \ n E^ be the spectral decomposition of M with all A„ two 
by two distinct. Now assume £ is an eigenvector of M with eigenvalue A„ ^ 0, 
satisfying E^( — (. Define rj by r\ r — ( r . Then one has 

(Mr)) r = M rsVs 

S 



= -(MtC) r 

= -XnVr- (67) 

Hence, n is an eigenvector of M with eigenvalue — A„. This implies that either 
A„ = or there exists m ^ n such that A m = — A„. In the latter case, m and 
n are matching indices and E^ projects on all vectors rj obtained by taking 
elementwise complex conjugation of all vectors in the range of 

Now choose an orthonormal basis , C'- 2 ) , • • • , m the range of £W and 

a corresponding basis Tf^\ rf 2 \ • • • , rf q ^ in the range of E^ n \ with n^P — ffl. 
Do this for all non- vanishing pairs of eigenvalues A m = — A„. Complement this 
with an orthonormal basis in the nullspace of M, if present. Collect all these 
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basis vectors as columns of a unitary matrix U. For a given £w) i n the range of 
£W is, with some abuse of notation, 

(U T MUSj) p = (U T MU) rj 

= A,<^)|C W ). (68) 

By construction, the latter vanishes for all but at most one value of r. This 
ends the proof. 

□ 

Finally, the above results are generalised to arbitrary square matrices. The 
argument is that found in the proof of [2], Lemma 1. 

Let be given a matrix M which is either symmetric or anti-symmetric. The 
matrix MM^ is hermitean and can be diagonalised by means of a unitary matrix 
U, i.e. WMM^U is diagonal. Let C = WM(W) T . Then C, like M, is either 
symmetric or anti-symmetric. In addition it satisfies (using that U'MM'U is 
diagonal and that M T = ±M) 

CC^ = U^MM^U = [U 1 MM^U) T = U t mUi{U t ^ = C^C. (69) 

This means that C is normal and that, by the previous propositions, there exists 
a unitary matrix V such that V T CV has on each row and each column at most 
one non- vanishing element. The proof of the two theorems then follows easily. 

Appendix B 

Here we present the calculation of the time average entanglement of the initial 
boson state 

V = |i,i,(V--,o> = &l4|o). (70) 

See Section gj 

The non-degenerate eigenvector of the two-particle Hamiltonian is 

4,m = _L 6 t (0 (o) )6 t W (o) ) | O) = _L_£ 6 t 6 t| 0) (71) 

It has eigenvalue 2. The projection of |1, 1, 0, 0, • • • ) onto this eigenvector is 
V^SJ^C 00 ) with = 2/N 2 . 

Introduce vectors £j , orthogonal to , determined by 



■0 ffl> +6 = ^ (O, +&- (72) 



Then one can write 



11,1,0,0,-.-) = l6t ( ^(o) )6 t (0 (o) ) | O) + 6 t (ei)6 t (6) |o) 
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+ -^t (0 (o) )6 t (6+6) | O ). 



(73) 

The projection of |1, 1, 0,0, ■ ■ ■ ) onto the (N — l) 2 -fold degenerate subspace 
equals & + (Ci)6 + (C 2 )|0). It is written as vV 11 ^ 1 ^ with 

P {11) = \\bHZi)b f (t,2)\0)\\ 2 

= (616) (616) + 1(6 16) I 2 

The projection of |1, 1,0, 0, ■ ■ ■ ) onto the remaining subspace equals 
-^&t (0 (o) )fet(a+6) |o). 



It is written as ^/p( 01 )?/;( 01 ) with 



P m - ^l|fc t (^ (o) )& t (6+6)|0)|| 2 
= ^II6+6II 2 

N \ N, 



Explicit expressions for the three eigenstates are 



^ = 



jk 
1 



1 



1 1 



(6 + 6) 



UJl ® 0J 2 + L0 2 ® Wl + ^<^ (0) ® <^ (0) 



-L</> (0) ® W 2 <g> (O) 



^(01) 



1 



1 



1 1 



v/iV 



(6+6) + (6+6)« 

(o) <g> (wi + w 2 ) + (wi + w 2 ) 



(74) 



(75) 



(76) 



(77) 



(78) 



The coefficients of the expansion of each of the vectors ^(°°) , i/>( n ) , 
into the basis vectors uij <g> cut can be written as 



A 



(oo) _ 1(1) 



'jk 



- 2 x jk > 
1 



1 



A (n) = 
jfe ^pW) 2NV2 

A (oi) = 1 1 



r (3) _ 2r (l) 



(79) 
and V (01) 

(80) 
(81) 

(82) 
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with 



2 

J 2 ) - 



x jk 

Vjk 



jl T Sj2)(8kl + Sk2) 

Sji + Ski + S j2 + 5 k 2, 
(Sji - 5 j2 )(5 k i - 5 k2 ) 



(83) 
(84) 

(85) 
(86) 



The matrices X^\ X^, X^ span a simple Jordan algebra of the spin factor 
type (see [HI], Section 2.9.7). The Jordan product is defined by 



One verifies that 



A*B = -(AB + BA). 



x m^ x (2) = 2X (2) 

X (3) *X^ 3) = 2X {3) + NX {2) + NX {1) 

X (D, X (2) = 2 (3) 
N 

X (1 KX ( V = 2X^+X^ 

x {2 Kx^ = X^ + 2X^ 

Y * Y = 2Y 

Y*X U) = 0, j = 1,2,3. 



(87) 



(89) 
(90) 

(91) 

(92) 
(93) 
(94) 
(95) 



There exists a representation of the Jordan algebra with the above product 
rules in R 2 + K + R, with the product rule 



(u, a, A) * (v, b, //) = (av + bu, (u\v) + ab, A/z). 

Let and be two unit vectors satisfying (u^\u^) ■ 
one can identify 



1,0) 



X& = ( U (2) ,l,0) 

X& = i(A>« W 2) ),4,0) 



Y = (0,0,2). 
With this representation is 



(96) 

-1 + 4/7V. Then 

(97) 
(98) 

(99) 
(100) 



v ^y A (oo) = 



2NV2 
1 

2Nt/2 
1 

NV2 



2X {1) - 2X {3) + NX {2) - NY 
(-(7V-2)u (1) ,iV-2,-2A r ) 

X (3) _ 2X d) 



(101) 



(102) 
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-±—UN-4)uM + Nu<®,Q,0). (103) 



It is now straightforward to calculate the squares 

P (00) (A (00) ) 2 = ^j(u (1) ,l,0) (104) 
1 

4AT2 



N 2 

p (n)( A (ii))2 = 1 (-(N -2) 2 u^\(N -2) 2 ,2N 2 ) (105) 



p (0i) (A (0i) )2 = J_ {0:N _ 2iQ y (106) 



Summing these relations gives 



a = r = l-(-(N-4)u^,N,2N). (107) 
Squaring again gives 

ct 2 =t 2 = — \-;(-2N(N - 4)u {1) ,N 2 + (N - 4) 2 ,4:N 2 ). (108) 
The trace of the matrix represented by (u, a, A) equals 2a + A. Hence, one finds 
Let us finally calculate A. Squaring again gives 

( p (oo )(A (oo ))2 ) 2 = JL( w Ml,0) (109) 

( p (ii) (A (n))2) 2 = _!_ ( _ (Ar _2)V 1 ),(iV-2) 4 ,2iV 4 ) (110) 

( p (01) (A (01) )2 ) 2 = i_ (0)(iV _ 2) 2 )0) . (m) 

This gives 

A = 1 - Tr (p(°°)(A( 00 )) 2 ) 2 - Tr (pC 11 ^ 11 )) 2 )' - Tr (p< 01 >(A< 01 >) 2 ) a 

— — — — 

2 + N~N 2 ~ + Jp~ly4' ( ' 
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